
;***************************************************************
;discription：计算像元某一个月的太阳总辐射量
;
;param1: nrows
;param2: ncols
;param3: dem_fid  dem图层
;param4: tmp_max 这个月每天最高气温，一共31天
;param4: tmp_min 这个月每天最低气温，一共31天
;param4: tmp_avr 这个月每天平均气温，一共31天
;param4: month 表示一年的月份
;param3: sol(output) 太阳总辐射量
;
;author: wangqiang( wensten )  2015-09-14  chengdu
;***************************************************************
PRO PRO_APAR_SOL, nrows, ncols, dem_fid, tmp_max, tmp_min, tmp_avr, month, sol=sol
  ;严格编译规则
  COMPILE_OPT IDL2
  ;二次开发模式初始化，恢复envi的sav文件
  ENVI,/RESTORE_BASE_SAVE_FILES
  ;初始化envi，保存日志文件
  ENVI_BATCH_INIT,log_file='sol.log'
  
  ;第一种，取定值：5852-6680mj/(m^2 * year) = 488-557mj/(m^2 * mon)
  ;西藏东部地区，取上限;注意上述单位是平方米，而下面计算是平方千米，需要进行单位转换
  ;sol=557 * 10^6
  
  ;第二种，动态计算：本文参考张佳华，等（2010，科学出版社）BEPS模型中总辐射计算公式
  num_days=31
  ;新建某一天的太阳辐射总量矩阵
  sol_day = FltArr(nrows,ncols,  1)
  PRO_APAR_SOL_DAY,nrows, ncols, dem_fid, tmp_max, tmp_min, tmp_avr, 1, sol_day=sol_day
;  FOR jday=(month-1)*num_days+1, month*num_days DO BEGIN
;    ;新建某一天的太阳辐射总量矩阵
;    sol_day = FltArr( ncols,nrows, 1)
;    PRO_APAR_SOL_DAY,nrows, ncols, dem_fid, tmp_max, tmp_min, tmp_avr, jday, sol_day=sol_day
;    ;月总量
;    sol=sol+sol_day
;    print,sol_day
;  ENDFOR
;  print,sol_day
  sol=sol_day
  envi_file_mng,id=dem_fid,/remove
  ;关闭envi二次开发模式
  ENVI_BATCH_EXIT
END



;***************************************************************
;discription：计算像元区域某一个天的太阳总辐射量
;
;param1: nrows
;param2: ncols
;param3: dem_fid 
;param4：tmp_max日最高气温
;param5：tmp_min日最低气温
;param6：tmp_avr日平均气温
;param7: day 表示一年中的第j天
;param8: sol_day(output) 日太阳总辐射量
;
;author: wangqiang( wensten )  2015-09-14  cdut
;***************************************************************
PRO PRO_APAR_SOL_DAY, nrows, ncols, dem_fid, tmp_max, tmp_min, tmp_avr, day, sol_day=sol_day
  ;严格编译规则
  COMPILE_OPT IDL2
  iproj = envi_get_projection(fid=dem_fid)
  oproj = envi_proj_create(/geographic)
  
  ;计算day在某个月中第几天
  indexofday = (day mod 31) - 1
  ;相对湿度
  humidity_cellvalue=0.6
  ;高程数据
  dims = [-1, 0, nrows - 1 , 0, ncols - 1]
  dem_arr=envi_get_data(fid=dem_fid, dims=dims, pos=0, xfactor=1.0, yfactor=1.0); get bands' value
  
  ;行方向(Y方向)，列方向(X方向)
  FOR yindex=0, nrows-1 DO BEGIN    
    FOR xindex=0, ncols-1 DO BEGIN
      ;获取点(xindex,yindex)处的高程值
      elev_cellvalue = dem_arr[yindex,xindex]
      IF (elev_cellvalue EQ -9999) THEN BEGIN
        sol_cellvalue=-9999
        sol_day[yindex,xindex]=sol_cellvalue  
      ENDIF ELSE BEGIN
        ;最高、最低、平均气温
        tmpmax_cellvalue=tmp_max[yindex,xindex]
        tmpmin_cellvalue=tmp_min[yindex,xindex]
        tmpavr_cellvalue=tmp_avr[yindex,xindex]
        IF (tmpmax_cellvalue EQ -9999) || (tmpmin_cellvalue EQ -9999)||(tmpavr_cellvalue EQ -9999) THEN BEGIN
          sol_cellvalue=-9999
          sol_day[yindex,xindex]=sol_cellvalue
        ENDIF ELSE BEGIN
          ;最高气温最低气温间隔
          tmp_range=tmpmax_cellvalue-tmpmin_cellvalue
  
          ;根据像素坐标(xindex,yindex)获取地图投影坐标
          envi_convert_file_coordinates, dem_fid, xindex, yindex, xmap, ymap, /to_map
          ;转换投影坐标到地理坐标
          envi_convert_projection_coordinates, xmap, ymap, iproj, longitude, latitude, oproj ;Lon = oxmap, Lat = oymap
          ;打印像素坐标、地图坐标、投影坐标和高程值
          ;print, xindex, yindex, xmap, ymap, longitude, latitude, elev_cellvalue
  
          ;大气顶层太阳辐射
          PRO_APAR_SOL_S0, latitude, day, s0=s0
          ;大气总透射率
          PRO_APAR_SOL_TCF,latitude, elev_cellvalue, tmpavr_cellvalue, humidity_cellvalue, tcf=tcf
          ;日长校正因子fd
          PRO_APAR_SOL_FD, latitude, day, fd=fd
          ;温度校正系数ft
          PRO_APAR_SOL_FT, tmp_range, tmpavr_cellvalue, ft=ft
          ;饱和水汽差校正系数fe
          PRO_APAR_SOL_FE, tmpmax_cellvalue, tmpmin_cellvalue, fe=fe
  
          ;计算像元的太阳辐射总量
          sol_cellvalue = tcf*fd*(1-ft*fe) * s0
          print, yindex, xindex, sol_cellvalue
          sol_day[yindex,xindex]=sol_cellvalue  
        ENDELSE
      ENDELSE
         
    ENDFOR
  ENDFOR
  
END



;***************************************************************
;discription：大气顶层太阳辐射
;
;param1: latitude 纬度,单位为角度
;param2: jday表示一年中的第j天
;param3: so(output) 大气顶层太阳总辐射量
;
;注意，在IDL中三角函数使用的弧度单位。角度转弧度：30*！DTOR ,弧度转角度：ACOS（0.5）*!DTOR
;author: wangqiang( wensten )  2015-10-3  chengdu
;***************************************************************
PRO PRO_APAR_SOL_S0, latitude, jday, s0=s0
  x=(2.0*!PI) * ( jday - 1.0/365 )
  ;太阳赤纬角，单位为角度
  delta= (180.0/!PI) * ( 0.006918-0.399912*cos(x)+0.070257*sin(x) - 0.006758*cos(2*x)+0.000907*sin(2*x) - 0.002697*cos(3*x)+0.00148*sin(3*x) )
  ;日地距离订正系数
  r=1.000110+0.034221*cos(x)+0.000719*cos(2*x)+0.000077*sin(2*x)
  ;太阳常数，单位： w/m^2
  I0=1367.0
  ;一天的长度，单位：小时h
  L=24
  ;日出日落时角，弧度表示
  omiga=acos( -tan(latitude*!DTOR) * tan(delta*!DTOR) )
  
  ;大气上界的太阳辐射. 注意,这里的单位是:W*M^-2(瓦/平方米)
  s0 = ( (L*I0)/( !PI* (r^2) ) )*( omiga*sin(latitude*!DTOR)*sin(delta*!DTOR) + cos(latitude*!DTOR)*cos(delta*!DTOR)*sin(omiga) )
  ;转换平方米到平方千米,注意：10^6=16960 , 10.0^6=1000000.0
  ;s0=s0*(10.0^6)
END


;***************************************************************
;discription：大气总透射率
;
;param1: latitude 纬度,单位为角度
;param2: elevtion 海拔
;param3: tmp 平均气温 
;param4: humidity 相对湿度
;param5: tcf(output) 大气总透射率
;
;author: wangqiang( wensten )  2015-10-3  chengdu
;***************************************************************
PRO PRO_APAR_SOL_TCF, latitude, elevtion, tmp, humidity, tcf=tcf
 
  ;与当地海拔有关的高程校正因子
  z = ( 1-( 2.2569/(10.0^5) ) * elevtion )^5.2553
  
  ;晴天条件下的透射率t0
  t0 =0.947 - ( (1.033/(10.0^5))*(latitude^2.22) )
  ;水汽透射率tv, 目前西藏研究区，空气干燥。注意：这个公式还有一个分支没有写完，它涉及到相对湿度。详细请见 《植被与生态遥感》 p154
  tv =0.9636 - (9.092/(10.0^5))*( (tmp+30.0)^1.8232 )
  ;大气分子和臭氧等引起的透射率ta
  ta =1.0
  
  ;大气总透射率
  tcf =( t0*ta*tv )^z
  
END


;***************************************************************
;discription：不同地区的日长校正因子fd
;
;param1: latitude 纬度，单位角度
;param4: jday 表示一年中的第j天
;param5: fd(output) 日长校正因子
;
;author: wangqiang( wensten )  2015-10-3  chengdu
;***************************************************************
PRO PRO_APAR_SOL_FD,  latitude, jday, fd=fd
  ;
  x=(2.0*!PI)*( jday - 1.0/365 )
  ;太阳赤纬角,单位为角度
  delta=(180.0/!PI) * ( 0.006918-0.399912*cos(x)+0.070257*sin(x) - 0.006758*cos(2*x)+0.000907*sin(2*x) - 0.002697*cos(3*x)+0.00148*sin(3*x) )
  ;日出日落时角,单位为弧度
  omiga= acos( -tan(latitude*!DTOR) * tan(delta*!DTOR) )  
  ;半天长
  H= 12.0 * (omiga/!PI)
  
  ;日长校正因子fd 
  fd= 1.0 / ( 1.0 - ( (H-!PI/4)^2 )/(2* H^2) )

END

;***************************************************************
;discription：温度校正系数ft
;
;param1: tmp_range 每日的最高和最低气温之间的温度范围
;param4: tmp 每日的平均气温
;param5: ft(output) 温度校正系数
;
;author: wangqiang( wensten )  2015-10-3  chengdu
;***************************************************************
PRO PRO_APAR_SOL_FT, tmp_range, tmp, ft=ft
  ;温度校正系数
  ft= 0.7 * 23.753 * tmp_range/( tmp+273.16)
END


;***************************************************************
;discription：饱和水汽差校正系数fe
;
;param1: tmp_max 每日最高气温
;param4: tmp_min 每日最低气温
;param5: fe(output) 饱和水汽差校正系数
;
;author: wangqiang( wensten )  2015-10-3  chengdu
;***************************************************************
PRO PRO_APAR_SOL_FE, tmp_max, tmp_min, fe=fe
  ;饱和水汽差校正系数
  fe= exp( 17.269*tmp_min/(273.3+tmp_min) ) / exp( 17.269*tmp_max/(273.3+tmp_max) )
END
